using Random, DelimitedFiles, PyPlot, StatsBase
Load reference data and define plotting (this will change after Przemek prepares better visualization)
loc = readdlm("loc.txt")
dist = readdlm("distkm.txt")
function plot_tour(t, loc)
t2 = [t; t[1]]
plot(loc[t2, 1], loc[t2, 2])
end
#This code is using route data files that need to be downloaded and unzipped from:
#https://github.com/pszufe/MapDatasets/raw/master/Ontario/walmart_routing.zip
pos = Dict{Int,Tuple{Float64,Float64}}()
for line in readlines("walmart_node_locations.txt")
els = split(line,"\t")
pos[parse(Int,els[1])] = (parse(Float64,els[2]),parse(Float64,els[3]))
end
routes = Matrix{Vector{Int}}(undef,size(dist)...)
fill!(routes,Int[])
f = open("walmart_routes.txt")
while !eof(f)
line = split(readline(f),"\t")
i,j = parse.(Int,line[2:3])
routes[i,j]= parse.(Int,split(line[12],"#"))
end
close(f)
k = first(keys(pos))
const MAP_BOUNDS = [[pos[k]...],[pos[k]...]]
for v in values(pos)
(v[1] < MAP_BOUNDS[1][1]) && (MAP_BOUNDS[1][1]=v[1])
(v[2] < MAP_BOUNDS[1][2]) && (MAP_BOUNDS[1][2]=v[2])
(v[1] > MAP_BOUNDS[2][1]) && (MAP_BOUNDS[2][1]=v[1])
(v[2] > MAP_BOUNDS[2][2]) && (MAP_BOUNDS[2][2]=v[2])
end
MAP_BOUNDS
#required installation for map vizualiztion
#using Conda
#Conda.runconda(`install folium -c conda-forge`)
using PyCall
flm = pyimport("folium")
function plot_tour2(t)
m = flm.Map()
for k=1:length(t)
i = t[k]
j = t[k<length(t) ? (k+1) : 1]
info = "Route $(k)-th <br>from=$i to=$j\n<br>distance=$(round(dist[i,j],digits=3))"
flm.PolyLine(
[pos[n] for n in routes[i,j]],
popup=info,
tooltip=info
).add_to(m)
end
for k=1:length(t)
i = t[k]
j = t[k<length(t) ? (k+1) : 1]
info = "Node $(k)-th <br>number=$i next=$j\n<br>\n$(round.(pos[routes[i,j][1]],digits=6))"
flm.Circle(
location=pos[routes[i,j][1]],
popup=info,
tooltip=info,
radius=2000,
color="crimson",
fill=true,
fill_color="crimson"
).add_to(m)
end
m.fit_bounds(Tuple.(MAP_BOUNDS))
m
end
Tour length calculation function common for all algorithms
function eval_tour(tour, d)
travel = d[tour[end], tour[1]]
for i in 2:length(tour)
travel += d[tour[i-1], tour[i]]
end
travel
end
Random initial tour generation common for all algorithms
init_tour(n) = randperm(n)
Function choosing one random neighbor used in local search and SANN
function pick_neighbor(tour, prob)
neighbor = copy(tour)
x, y = rand(eachindex(tour), 2)
if x != y
if rand() < prob
x,y = minmax(x,y)
reverse!(@view neighbor[x:y])
else
tmp = neighbor[x]
if x < y
neighbor[x:y-1] = @view neighbor[x+1:y]
else
neighbor[y+1:x] = @view neighbor[y:x-1]
end
neighbor[y] = tmp
end
end
return neighbor
end
Local search routine
function localsearch(dist, iter, prob)
t = init_tour(size(dist, 1))
v = eval_tour(t, dist)
best_t, best_v = t, v
for i in 1:iter
newt = pick_neighbor(t, prob)
newv = eval_tour(newt, dist)
if newv < v
t, v = newt, newv
if v < best_v
best_v, best_t = v, t
end
end
end
(best_v, best_t)
end
alg_localsearch(;iter=10^7, prob=0.8) = (x -> localsearch(x, iter, prob), "Local search (iter=$iter, prob=$prob)")
Simulated annealing routine
function sann(dist, T0, Tf, α, iter, prob)
T = T0
t = init_tour(size(dist, 1))
v = eval_tour(t, dist)
best_t, best_v = t, v
while T >= Tf
for i in 1:iter
newt = pick_neighbor(t, prob)
newv = eval_tour(newt, dist)
if newv < v || rand() < exp(-(newv - v) / T)
t, v = newt, newv
if v < best_v
best_v, best_t = v, t
end
end
end
T *= α
end
(best_v, best_t)
end
alg_sann(; T0=1, Tf=0.001, α=0.9, iter=round(10^7 * log(α)/ log(Tf/T0)), prob=0.8) =
(x -> sann(x, T0, Tf, α, iter, prob), "SANN (T0=$T0, Tf=$Tf, α=$α, iter=$iter, prob=$prob)")
Genetic algorithm routine
function init_tour_ga(d)
t = randperm(size(d, 1))
(eval_tour(t, d), t)
end
function mutate(tour, d)
x, y = minmax(rand(eachindex(tour[2]), 2)...)
reverse!(@view tour[2][x:y])
return (eval_tour(tour[2], d), tour[2])
end
function crossover(t1, t2, d)
x, y = minmax(rand(eachindex(t1[2]), 2)...)
child = similar(t1[2])
used = Set{Int}()
for i in x:y
child[i] = t1[2][i]
push!(used, t1[2][i])
end
pos = 1
for j in eachindex(t2[2])
if pos == x
pos = y + 1
end
if !(t2[2][j] in used)
child[pos] = t2[2][j]
pos += 1
end
end
(eval_tour(child, d), child)
end
function ga(dist, popsize, iter, mutprob)
w_replace = Weights(0:popsize-1)
w_parent = Weights(popsize:-1:1)
pop = [init_tour_ga(dist) for i in 1:popsize]
sort!(pop)
best_v, best_t = pop[1]
for i in 1:iter
replace_idx = sample(1:popsize, w_replace)
pop[replace_idx] = crossover(pop[sample(1:popsize, w_parent)],
pop[sample(1:popsize, w_parent)],
dist)
for j in 2:popsize
if rand() < mutprob
pop[j] = mutate(pop[j], dist)
end
end
sort!(pop)
if pop[1][1] < best_v
best_v, best_t = pop[1]
end
end
(best_v, best_t)
end
alg_ga(; popsize=120, iter=3*10^5, mutprob=0.01) =
(x -> ga(x, popsize, iter, mutprob), "GA popsize=$popsize, iter=$iter, mutprob=$mutprob")
Tabu search routine
function list_neighbors(tour, d)
neighbors = Vector{typeof(tour)}()
n = length(tour)
for x in 1:n-1, y in x+1:n
neighbor = copy(tour)
reverse!(@view neighbor[x:y])
push!(neighbors, neighbor)
end
sort!([(eval_tour(nei, d), nei) for nei in neighbors])
end
function tabusearch(d, tabulength, iter)
t = init_tour(size(dist, 1))
v = eval_tour(t, d)
best_t, best_v = t, v
tabulist = [Int[] for i in 1:tabulength]
tabupos = mod1(2, tabulength)
tabulist[tabupos] = t
for i in 1:iter
neighbors = list_neighbors(t, d)
found = false
for (neiv, nei) in neighbors
if !(nei in tabulist)
found = true
tabulist[tabupos] = nei
tabupos = mod1(tabupos+1, tabulength)
t, v = nei, neiv
if v < best_v
best_v, best_t = v, t
end
break
end
end
if !found
break
end
end
(best_v, best_t)
end
alg_tabusearch(; tabulength=100, iter=1000) =
(x -> tabusearch(x, tabulength, iter), "Tabu search tabulength=$tabulength iter=$iter")
Algorithm comparison
function run(;loc, dist, alg)
@time res = alg[1](dist)
#plot_tour(res[2], loc)
#title("$(alg[2]); cost: $(res[1])")
plot_tour2(res[2])
end
Random.seed!(2019);
run(loc=loc, dist=dist, alg=alg_localsearch())
run(loc=loc, dist=dist, alg=alg_sann())
run(loc=loc, dist=dist, alg=alg_ga())
run(loc=loc, dist=dist, alg=alg_tabusearch())